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Abstract 

We consider ensembles of sine-coupled phase oscillators consisting of subpopulations 
of identical units, with a general heterogeneous coupling between subpopulations. 
Using the Watanabe-Strogatz ansatz we reduce the dynamics of the ensemble to a 
relatively small number of dynamical variables plus microscopic constants of motion. 
This reduction is independent of the sizes of subpopulations and remains valid in 
the thermodynamic limits, where these sizes or/and the number of subpopulations 
are infinite. We demonstrate that the approach to the dynamics of such systems, 
recently proposed by Ott and Antonsen, corresponds to a particular choice of micro- 
scopic constants of motion. The theory is applied to the standard Kuramoto model 
and to the description of two interacting subpopulations, exhibiting a chimera state. 
Furthermore, we analyze the dynamics of the extension of the Kuramoto model for 
the case of nonlinear coupling and demonstrate the multistability of synchronous 
states. 

Key words: Coupled oscillators, oscillator ensembles, Kuramoto model, nonlinear 
coupling 
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1 Introduction 



A model of all-to-all, or globally coupled limit cycle oscillators explains many 
natural phenomena in various branches of science. The applications range from 
the description of the collective dynamics of Josephson junctions PQ, lasers [2J, 
and electrochemical oscillators [3 J to that of pedestrians on footbridges [4"f5] . 
applauding persons in a large audience [6], cells, exhibiting glycolitic oscil- 
lations [7,8,9] , neuronal populations [ID] , etc. Externally forced or feedback 
controlled globally coupled ensemble or several interacting ensembles serve as 
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models of circadian rhythms, normal and pathological brain activity, interac- 
tion of different brain regions, and many other problems |Hfl2fl3fl4fl5fl6fl7] . 
Many aspects of the ensemble dynamics, especially those related to inhomo- 
geneity [T5fl9f2U] or nonlinearity of coupling [2"T|2"2"] , temporal dynamics of the 
collective mode [23~f24"] , different frequency distributions [25|26] , and clustering 
[27|28|29] remain in the focus of the current research activity. 

Ensembles of weakly interacting units are successfully treated within the 
framework of phase approximation [3T)|31|I32|I33|I34] . Most popular is the Ku- 
ramoto model of sine-coupled phase oscillators. This model explains self- 
synchronization and appearance of a collective mode (mean field) in an ensem- 
ble of generally non-identical elements; the transition to synchrony occurs at 
a certain critical value of the coupling constant that is roughly proportional to 
the width of the distribution of natural frequencies [3T)f3~T] . With the further 
increase of coupling, more and more oscillators join synchronous cluster, so 
that the amplitude of the mean field grows as a square root of supercriticality. 
It is instructive to interprete this transition as follows: the non-zero mean field 
forces individual units and entrains at least a part of them; these entrained 
units become coherent, thus yielding a non-zero mean field. A quantitative con- 
sideration, based on this argument and first performed by Kuramoto [3U|3T] . 
provides the amplitude and frequency of the stationary solution. References 
to many further aspects of the Kuramoto model can be found in [35.36.37j. 

An extension of the Kuramoto model for the case of nonlinear coupling has 
been suggested in our recent publications [2"Tf2"2"] . where the most simple case 
of identical units has been treated. Nonlinearity in this context means that the 
effect of the collective mode on an individual unit depends on the amplitude 
of this forcing, so that, e.g., the interaction of the field and a unit can be 
attractive for weak forcing and repulsive for strong one. This can lead to 
nontrivial effects like a destruction of a completely synchronous state and 
appearance of partial synchrony in an ensemble of identical units. Moreover, 
in this state the frequencies of the collective mode and of oscillators can be 
different and incommensurate. The analysis of the nonlinear model for the 
case of an ensemble with a frequency distribution is still lacking and will be 
performed below. This analysis is based on the extension of the Watanabe- 
Strogatz (WS) theory [3"51l3"9"] to the case of nonidentical oscillators, suggested 
in our brief communication |40j . 

For a population of identical oscillators, a full dynamical description of the Ku- 
ramoto model can be obtained with the help of the WS ansatz which reduces 
the dynamics to that of three macroscopic variables plus constants of motion. 
This works even for a nonlinear coupling [22] . The main idea of this paper is in 
extending the WS ansatz on an ensemble of nonidentical oscillators, treating it 
as a system of coupled subpopulations, each consisting of identical oscillators. 
Each subpopulation can be then described by three WS variables, whereas 
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the full system is described by a system of coupled WS equations. A descrip- 
tion of an ensemble with a continuous frequency distribution is then obtained 
by performing a thermodynamic limit. The extended WS approach is applied 
to a population of oscillators with a Lorentzian frequency distribution (the 
standard Kuramoto model), and to a population of identical oscillators with 
inhomogeneous coupling [18]. In particular, within this framework we estab- 
lish a relation between the WS approach and the recent Ott-Antonsen (OA) 
ansatz [23,24], which yields, under certain assumptions, a dynamical equation 
for the evolution of the mean field. In this paper we discuss in detail how 
the WS theory can be applied to populations of non-identical units and apply 
this approach to describe the dynamics of oscillator ensembles with global 
nonlinear coupling, for the case of Lorentzian distribution of frequencies. 

The paper is organized as follows. We discuss the main model in Section |5J 
Extension of the WS theory to the case of nonidentical oscillators is presented 
in Section |3j In Section [4] we discuss a relation between the WS theory and 
the OA ansatz (cf. a recent paper by Marvel, Mirollo, and Strogatz |JT] for a 
related discussion) . This theory is applied to describe the dynamics of linearly 
(Section [5]) and nonlinearly (Section E]) coupled ensembles; here we also sup- 
port the theory by numerics. In particular, we demonstrate the differences of 
the dynamics of the full equation system and that confined to the OA reduced 
manifold. We summarize and discuss our results in Section [7J 



2 Hierarchically organized population of oscillators 

In this Section we introduce a model of hierarchically organized populations 
of oscillators and describe it in terms of microscopic equations of motion. The 
main idea is to treat an ensemble of nonidentical oscillators as a collection of 
subpopulations of identical oscillators. 

We start by introducing an ensemble of nonidentical phase oscillators, charac- 
terized by phase variables (fik and generally time-dependent frequencies 0Jk{t), 
where k = 1, . . . , N is the oscillator index and N is the ensemble size. Each 
oscillator generally interacts with all other units and is subject to external 
fields. The effect of these interactions can be represented as some effective 
forcing and therefore each unit is described as a driven oscillator 



where variables and characterize the amplitude and phase of the force. 
It is convenient to introduce a complex force Hk = Af.e^ h and re-write Eq. (TjQ) 
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as 

^=u k + lm(H k e-^) . (2) 

In many particular cases considered below the force H is calculated in some 
mean-field manner from the state of the whole population or some subpopula- 
tion(s). However, for a while we prefer to consider H, as well as u, as arbitrary 
functions of time; for brevity we skip this dependence in the notations. Note 
that generally H can include random component (s) which describe a forcing 
by common noise. We emphasize that Eqs. (j2J) do not represent the most 
general model of coupled phase oscillators, because the high order harmonic 
components of the phase ~ e~ tn ^ k , w jth n > 1, do not appear on the right 
hand side of Eqs. (J2J). 



external 




Fig. 1. Hierarchically organized ensemble (left panel) can be considered as a collec- 
tion of subpopulations labeled by index a = 1, . . . , M, with a generally bidirectional 
asymmetrical coupling between subpopulations. Generally, the subpopulations are 
exposed to external field(s). Each subpopulation (right panel) of the oscillator en- 
semble consists of identical oscillators driven by a common force. This force results 
from the action of all other subpopulations and external fields and from the cou- 
pling within the subpopulation, which is usually taken to be of a mean field type. 
This means that outputs of all oscillators create the mean field that acts back on 
these oscillators. The crucial point is that the forcing is formed either via linear or 
nonlinear transformation of acting mean fields. 

In general, all oscillators described by fl2]) have different dynamics, and such 
a system cannot be further reduced. Such a simplification is possible for an 
ensemble of identical oscillators, where the WS theory holds which states that 
dynamics of a population of identical oscillators of type (J2J) subject to a com- 
mon forcing is effectively three-dimensional (this is discussed in details in 
the next Section). In order to make use of the WS ansatz to an ensemble of 
nonidentical oscillators, we assume that this ensemble consists of groups (sub- 
populations) of identical units (Fig. [fl cf. [42j ) . Accordingly, we re-label the 
oscillators grouping them into M subpopulations, denoted by index a. Each 
population contains N a identical units with frequency uj a which are driven by 
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a common force H a . The equations now read: 



(a) 




k 



u a + Im H a e 



(3) 



dt 



where k — 1, . . . , N a and a — 1, . . . , M, with an obvious relation J2 a =i N a = N. 
Note that two subpopulations can have the same frequencies u a = Ub but differ 
by the force H a ^ Hb, or they can have different frequencies but the same force, 
or they can differ both by their frequencies and the forces. 

Certainly, such grouping of oscillators into subpopulations is not always pos- 
sible (if we exclude a trivial case when each subpopulation contains one unit 
and M = N). However, in many particular cases M <C N and the description 
of the dynamics simplifies essentially: as we will see in the next Section, each 
group is described by three WS variables and instead of N equations (j3J) we 
have to deal with M coupled systems of three WS equations each. Note that 
WS ansatz is also valid for infinitely large populations of identical oscillators. 
Thus, if some or all of subpopulation sizes N a — > oo, we still obtain a 3M- 
dimensional description of the collective dynamics. As discussed in the next 
Section, the idea to consider an ensemble as a collection of subpopulations is 
also fruitful for treating a thermodynamic limit iV — » oo with a continuous 
distribution of oscillator frequencies. 



3 Oscillator populations in external fields 

In this Section we discuss the dynamics of large populations of oscillators 
subject to an arbitrary external force. First, we briefly present the WS theory 
which treats ensembles of identical oscillators. We re- write the WS equations in 
new notations, what makes them more convenient for the consequent analysis. 
Next, we extend the WS theory to the case of a finite number of interacting 
groups. Finally, performing a thermodynamic limit, we obtain a description 
of an ensemble with a continuous frequency distribution. 



3.1 Identical oscillators 



First we discuss in details dynamics of a population of identical oscillators 
subject to an arbitrary forcing, which is, however, common for all oscillators. 
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3.1.1 WS equations 



The main result of the seminal papers by Watanabe and Strogatz [38|l39j is 
that a population of iV > 3 sine-coupled phase oscillators with any time- 
dependent common frequency u(t), driven by an arbitrary but common force 
Hit), 

^ = u(t) + Im (H(t)e-^) , k = 1, . . . , N , (4) 

can be completely described by three global variables p(t), $(t),and plus 
constants of motion ip^ that are determined by initial conditionso The original 
phases can be recovered by means of the time- dependent transformation 

p *<t> k _ p i$ tL 1 (K) 



p e iWk-*) + 1 



or, equivalently, 



1-p ( iP k -* \ 
tan {-^) = —p tan {-^) ' (6) 

see Fig. [2] for illustration^] 

Here ipk are the constants of motion. Notice that only N — 3 of them are 
independent; see the discussion below. The time-dependent functions p(t), 
and are the global amplitude and phase variables, respectively. In 
the following we denote them as WS variables, their meaning is discussed 
later. For the transformation (JSJ) to be consistent with equations of motion 
(jH), these variables have to obey the WS equations 



dp 1 — p 



2 



Re(He- 1 *) , (7) 



dt 2 
_ = w + _l m( ff e ; 



d^ 1 - p 



2 



dt 2p 



Im(He-**) . (9) 



We emphasize that we use the global variables slightly different from those 
originally used by WS [39J and write the equations in a different form. The 
relation to the original form can be found in Appendix |XJ there we also demon- 
strate the equivalence of the transformations (jSJ) and (jSJ). 

For a further analysis it is convenient to introduce a combination of two WS 
variables z = pe 4 *. We call this complex variable the bunch amplitude. Intro- 



1 There is a restriction: an initial state of the ensemble cannot contain too large 
clusters, see [39J for details. 

2 Note that transformation ip — > <j) has a form of the Mobius transformation |41j . 
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Ipk-^ 



Fig. 2. Illustration of the transformation ipk fat see Eqs. (|5l6p . here for p = 0.5 
(cf. Fig. 3 in [39]). 



ducing also the phase shift a = $ — we write the WS equations ([319]) in an 
equivalent form] 



dz 
~dt 
da 
~dt 



1 z 2 
iujz + -H H* 

2 2 



w + lm(z*H) . 



(10) 
fill 



Constants of motion 

Transformation (jSJ) from original variables 0^ to WS variables p, $, ^ (or 2; 
and a) and ^ yields an overdetermined system. Hence, we have to impose 
3 constraints on the constants ipk] this is discussed in detail in [39]. It is 
convenient to choose two constraints as follows: 

N 

e^ fc = . (12) 

k=i 

The third constraint on is somehow arbitrary, it only fixes the common 
shift of ipk with respect to ^. It can be taken, e.g., as J^k^k — (this implies 
that constants are defined in the [— 71, 71] interval) [39]. Another convenient 
choice is cos(2ipk) = 0. The imposed constraints allow one to determine 
unambiguously the new variables p(0), ^(0), $(0) and constants ipk from the 
initial conditions </>fc(0) and vice versa, this is discussed in details in Ref. [3~9~| . 



3 We note here that Eq. fllOf) coincides with the Ott-Antonsen equation [23] for 
the dynamics of the complex mean field. However, in the Ott-Antonsen ansatz it 
appears without Eq. (jllh . The relation between the OA ansatz and the WS theory 
is treated in details in Section |H 



3.1.3 Meaning of the WS variables 

In order to discuss the physical meaning of the WS variables we compare them 
with the complex mean field, or the Kuramoto order parameter 

N 

Z = N" 1 Y, = re i@ , (13) 

k=l 

where r and G are the amplitude and phase of the mean field, respectively. [The 
comparison here is qualitative; a quantitative relation is given in Section HI] 

The WS amplitude variable p is roughly proportional to the mean field am- 
plitude r. Indeed, if p = 0, then from Eq. (jSJ) with account of Eq. (TT2|) . we 
obtain r = 0. Similarly, Eq. fl5]) shows that for p = 1 all <pk are equal, what 
yields r = 1. For intermediate values < p < 1 the relation between p and r 
generally depends also on \1/ and ipk- 

The WS phase variable $ characterizes the position of the maximum in the 
distribution of phases and is close to the phase of the mean field 0. They 
coincide for p — r — 1; for < p < 1, $ is shifted with respect to 9 by a 
factor, dependent on p, \F 




Fig. 3. Illustration of the meaning of the WS variables. The global amplitude p is 
related to the width A of the distribution of phases; p = if this distribution is 
uniform and p = 1 if it collapses to a 5-distribution. Thus, p is roughly proportional 
to the mean field amplitude r. Angle variable describes the position of the hump 
in the distribution; therefore it roughly corresponds to the phase G of the mean field. 
Angle variable Vl/ characterizes the motion of individual oscillators with respect to 
the hump. 

The second WS phase variable \1/ determines the shift of individual oscillators 
with respect to $. From Eqs. ( |8ll9l) we obtain a useful relation: 

* = 1^4(6 (14) 
1 + p l 
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As can be seen from Fig. [2] and Eq. ([6]), <pk decreases by 2n if \1/ grows by 2n. 
Hence, the oscillator frequency is 

COosc = (k) = ($) " W • (15) 



3.1.4 Dynamics of WS variables 

Here we discuss possible solutions of Eqs. (EH]) for given ou(t), H(t). First, we 
note that these equations represent a skew system: the variable \1/ does not 
enter Eqs. f JTfHj) . Hence, one has to solve the two-dimensional system f!7|8|) and 
then obtain \l/(t) via Eq. (TH|) . Next, the system (J7J|9]) obviously admits a fully 
synchronous solution p = 1, describing a synchronized cluster of identical 
oscillators: all in the same state. As it follows from Eq. in this case \l/ 
is an arbitrary constant and the global dynamics is described solely by an 
equation for $. However, we cannot analyze stability of this solution for general 
functions u(t), H(t). We mention here a seemingly trivial state when Hq = 0. 
Then p is arbitrary, $ = oj, and \l/ = const. This solution appears in the 
nonlinear model treated below in Section EJ 

It is instructive to find solutions of Eqs. ([312]) for an important particular case 
u = const and H(t) = H e tut . For this case of a harmonic forcing, Eqs. (JTHH]) 
take the form: 

^ = i^#ocos(i/*-$), (16) 

^ = w + i±^flosin(^-$), (17) 
at 2p 

aim 1 - p 2 

— = -^-H Bm{vt - *) . (18) 

differences 
this system as an autonomous one 



7T 

Introducing the phase differences A = — + $ — vt and a = $ — ^ we rewrite 



f = ^osmA, (19) 

^ = w _ z/ + 1±Z/7 oCOS a (20) 
at 2p 

— = ou — v — H p cos A . (21) 
at 



Remarkably, the system (I19|20l) is reversible, as it remains invariant under a 
transformation A — > —A, t — > —t. The system possesses two types of solutions, 
in dependence on the detuning \cu — v\. For \u — v\ < H system ( IT91I201) has one 
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attractive and one repelling steady state solution: these fixed points are located 
symmetrically with respect to the involution A — > —A and have coordinates 
p = 1 and A = ± arccos , respectively. For the case \u — v\ > H the 
system has no attractors, and the exact solution, described in Appendix [B], is 
represented by a family of closed orbits around a marginally stable equilibrium 
with coordinates A = and po (see Eq. fIB.lj) ). For this equilibrium solution, 
\l/ rotates with a constant frequency, given by Eq. flB.2j) . Except for this special 
case, the variables p and A also oscillate in time; moreover, the variable \l/ 
possesses an additional frequency so that the full dynamics is quasiperiodic 
(see Appendix [Bl. 

In summary, the behavior of the ensemble of non-interacting oscillators in a 
common field "follows" the dynamics of an individual oscillator. If the latter 
is entrained by the force then the ensemble is also fully synchronized; this 
corresponds to an attractive fixed point solution of system (I19ll20p . If an indi- 
vidual oscillator is not locked to the field, then there is no effective dissipation 
in the ensemble and system (I19|I20I) has no attractive solutions; as a result 
the collective variables oscillate. Only for a special preparation of initial con- 
ditions these variables are p = po and A = 0, i.e. they do not vary in time. 
This dynamics of the collective variables has a clear physical meaning if we 
interprete it as a dynamics of macroscopic characteristics of the population 
distribution. If the oscillators are locked, the population distribution collapses 
to a 5-function, what corresponds to p = 1, and the only relevant quantity is 
the position of the cluster $. If the oscillators are not locked, then each of them 
rotates non-uniformly with respect to the external field, and the distribution is 
generally "breathing" - unless the initial conditions correspond to the steady 
distribution, i.e. to the marginal equilibrium point p = p and A = 0. From 
this picture follows an important fact: the time averages along all possible tra- 
jectories of ( TT9ll2Tj) are equal. Indeed, average of a macroscopic variable can 
be obtained by averaging over time average for individual oscillators, and the 
latter quantities are the same because the oscillators are identical. 



3.2 Ensemble of nonidentical units as a hierarchical population 



Consider now an ensemble which can be viewed at as a collection of M subpop- 
ulations of identical units. Let index a = 1, . . . , M label the subpopulations 
and let each subpopulation consist of N a > 3 oscillators. Again, we assume 
that the force H a acting on a subpopulation a equally effects all oscillators 
of this subpopulation. Hence, the dynamics of each subpopulation can be de- 
scribed by means of the WS ansatz. For the a-th subpopulation we write, 
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similar to fllOjlip 



iu a z a + l -H a - ^H* a , (22) 
co a + lm{z* a H a ) . (23) 

For each subpopulation we have to specify constants of motion tp a ,k, k = 
1, . . . , N a ; each set of these constant obeying the same three additional con- 
straints, see Eq. fTl2|) and the paragraph following it. 

Hence, the full hierarchical system is described by 2M equations for M com- 
plex variables z a and M real variables a a , plus N — 3M constants of motion. 
(Certainly, we can use the WS equations in the real form to obtain 3M 
equations for 3 real variables.) 



dZg 

dt 
da a 



3.2.1 Treating individual oscillators and clusters 

Having introduced the coupled WS equations for several interacting subpop- 
ulations we come back to the limitation of the WS theory: the number of 
oscillators should be larger than three, and an initial configuration of a sub- 
population cannot have too large clusters of fully identical oscillators. To over- 
come this, it is sufficient to notice that Eq. (T5]) for an individual oscillator can 
be also written in form ([7J|9]) if we set p = 1, $ = 0, and \1/ as an arbitrary 
constant. In this case in (f22|) we have \z\ = 1 and the two equations of system 
( 122|23|) are in fact equivalent. Thus, individual oscillators not belonging to 
large groups can be also described by Eqs. f)22|23p . 

The same idea helps to treat large clusters in subpopulations. Suppose that 
a subpopulation a with N a elements does have a fully synchronized cluster 
formed by Nff> > N a /3 oscillators. To treat such a group, we split this sub- 
population into two, labeled by a' and a", respectively. (It means that we have 
now M + 1 subpopulations.) The first one is formed by Nff' elements of the 
cluster, the second one is formed by other iV a — iV^ c ) oscillators. The first sub- 
population can be treated like one oscillator: it has p a > = 1 and is described by 
one equation for the phase $ a / (see Section [3. 1.41) : the second one is described 
by three WS equations, so that altogether we have four equations for this 
subpopulation. If it turns out that the population a" also contains a majority 
cluster, then again the elements of this cluster can be considered as a separate 
subpopulation, and so on. 
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3.3 Infinitely large populations 



Most applications of the theory treat infinitely large ensembles. Here we con- 
sider separately two cases: (i) the number of subpopulations remains finite, but 
the subpopulations are infinitely large and (ii) the number of subpopulations 
is infinite. 

3.3.1 Thermodynamic limit I: Finite number of subpopulations 

If all or some of M subpopulations are infinitely large, they still can be de- 
scribed by Eqs. ( 12211231) . However, any infinitely large subpopulation is now 
described not by a finite set of constants of motion ipa,k, but by the distribu- 
tion functions cr a (if)) (see [39]). Equations ( !T2|) take the form 

f a a (^)e^# = . (24) 

J — 7T 

The third constraint for the constants of motion is also written via an integral, 
e.g., J^ n a(ip) cos 2ipdip = 0. Correspondingly, the phases are described by a 
distribution density w a ((f)). 

3.3.2 Thermodynamical limit II: Infinite number of subpopulations 

Now we assume that the number of subpopulations M — > oo. Hence, we 
substitute the subpopulation index a by a continuous variable, say x. In most 
applications this variable will be associated with the frequency, and in this 
way we will be able to describe ensembles with a continuous distribution of 
frequencies, however at the moment we prefer not to specify the meaning of 
x. (Generally, x can be also a vector variable.) 

Thus, the WS variables z, a or p, $, as well as the forcing H become 
functions of x and t and the WS equation system takes the form: 

1 z 2 
iu(x,t)z + -H(x,t)-—H*(x,t), (25) 

u(x,t)+lm(z*H(x,t)) . (26) 

In performing this limit we can consider the subpopulations as finite or infinite. 

We note that generally, the description of the system in this limit remains 
infinitely dimensional. However, it is simpler than the original one because for 
each value of the continuous variable x, which can have, e.g., the meaning of 
frequency, we have only three WS variables. Moreover, as is discussed below, at 



dz(x, t) 

dt 
da(x, t) 

dt 
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least for the case of global coupling with the Lorentzian frequency distribution, 
the description becomes really low- dimensional. 



3.3.3 Direct WS reduction for a system with continuous distribution of pa- 
rameters 



Above we have treated an ensemble with a continuous distribution of param- 
eters as a thermodynamic limit of a hierarchically organized populations. As 
the number of elements at each value of the continuous parameter is arbitrary, 
it appears instructive to apply the WS ansatz directly to an equation for the 
distribution density. Our derivation, presented below, is heavily based on the 
derivation given by WS (see Ref. [2S]), where they treated the case of identical 
oscillators. 



We start with the continuity equation which expresses the conservation of the 
number of oscillators: 

dw d 

«- + 5f <«.) = <>, (27) 

where w(x,<fi,t) is the distribution density function. The velocity v = <p is 
determined by the microscopic equation of motion 

v = u(x,t) + lm \H(x,t)e^ 

Following the idea of Watanabe and Strogatz |39j we demonstrate that, with 
the transformation to the WS variables p, $, \1/ and ip, the time-dependent 
density w(x,<j>,t) is transformed into a stationary density a(x,ip). 

We perform the following variable substitution in the continuity equation: 



t, </>, x —> r = t, ip = if>(x, <f>,t), y = x . 



The relation between the densities in old and new variables takes the form: 

d(y ip t) dip 
w(x, <p, t) = a(y, ip, t) - - = a(x, ip,r)— . (28) 

0{X, (f), t) 0(f> 

Writing the continuity equation in new coordinates (see Appendix O . we 
obtain: 

dw d . da dip da j dip ( dip dip \ 1 

= ~dt + W WV ' = !h~£kp + d^\d^\dt +V d^)j 

2 . (29) 
d ( dip \ d { dip\ ( dip dip\ ( dip\ dv 

'° < Ih { d^p 1 + !hp \ ~d~4> ) V ~dt + V ~d4> ) + [ity I Ihp 
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It AppendixOwe show that both expressions in curly brackets vanish provided 
$(x, t), ty(x,t), p(x,t) obey : 



dp(x, t) 1 — p 2 
~~dt ~ 2 



Re{H(x,t)e-^) , 



d$(x,t) 
dt 



uj(x,t) + 




Im(#(M)e~'*) , 



(30) 



2p 



dV(x,t) _ 1-p 2 
dt ~~ 2p 



lm{H{x,t)e-^) . 



da 



This implies that — = and, thus, a(x,ip) is a stationary distribution. 

Hence, the transformation to WS variables indeed results in a low-dimensional 
description of the dynamics (three global variables and function a for each 
value of x). Equations (|30|) are equivalent to Eqs. (I25ll26p . 



4 Linking the Watanabe-Strogatz theory and the Ott-Antonsen 



In this Section we relate WS variables to the complex mean field (which is the 
Kuramoto order parameter, see Eq. (fl3|) ). and to the generalized Daido order 
parameters. We demonstrate that particular solutions of the WS equations for 
the uniform distribution of constants of motion are equivalent to the solutions 
obtained with the help of the Ott and Antonsen ansatz [2321], see also |4"3"f4"l] 
Next, we discuss properties of the OA solution for singular and continuous 
distributions of oscillator frequencies. Note, that a relation between the WS 
and OA theories has been also recently discussed by Marvel, Mirollo and 
Strogatz HI]. 

4-1 WS variables versus order parameters 

For a hierarchically organized population we can define the mean field for each 
subpopulation; we call this quantity local mean field. If the subpopulation 
is finite, its mean field Z a is computed via Eq. ffl3|) . For an infinitely large 
subpopulation it can be computed as 



where w a ((p) is the probability distribution density of phases in the subpopu- 
lation a (cf. Section [3.3. 1 j) . 



ansatz 




(31) 
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Let us relate Z a and WS variables. [For simplicity of presentation we omit the 
index a in the following equations in this subsection.] With the help of Eq. ((SJ) 
we obtain: 



N 

Z = re i0 = N- 1 2 e irpk = pe i<& 7 (p, *) = zj{z, a) , (32) 
fc=i 

where 

^ 1 _|_ 

7fa = 1+ % (ltt -„ , PS) 

or, using another set of variables, 

7(3,01 ) = iV- 1 V^-U 77 — . 34 

If the subpopulation is infinite, then the summation is substituted by the 
integration and we have, e.g., instead of (1331) 

r-K i _i_ /i- 1 ( =*(V'-'i') 
= 1+ P pe „,-„ (35) 

We see that in general a relation between the WS variables and the order 
parameter is rather complex and contains not only the macroscopic variables 
p, $ but also all microscopic constants ipk (or, for an infinite subpopulation, 
depends heavily on the distribution 

Now we discuss an important particular case when the order parameter Z can 
be expressed only through two WS variables p and $. For this purpose we 
re-write the function 7 (see Eq. 1331) as a series: 



N 

7 (p,¥) 

fc=i 



00 , 



(i + p ~VW' fc -*))^(_ pe ^ fe -* 

f^-^e - **) ("at- 1 53 e ^ fc ) + 
\z=o / V fe=i / 

co / N 



Z=0 \ fc=l 

00 , 00 



«=o z=o 

where the coefficients 



A? 



d = N- 1 e^ fc (36) 

k=l 

are the amplitudes of the Fourier harmonics of the distribution of constants 
of motion ip^. Using that C\ — due to Eq. ( fl2l) or Eq. ( I24"l) . respectively, we 
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finally obtain 

oo 

7 = l + (l-p- 2 )£<3(-pe- i *)' • (37) 

1=2 

The same series representation for 7 is obtained in the thermodynamic limit 
when 7 is computed via integration (see Eq. fl35|) ). In this case the coefficients 
are computed according to 

J —TV 



The crucial observation is that Eq. (13TI) essentially simplifies and we obtain 
simply 7 = 1 if all the amplitudes of the Fourier harmonics vanish, i.e. if Ci = 
0. For an infinitely large system this is obviously true for a uniform distribution 
of constants of motion, i.e. for o'(u) = 2n~ 1 . For a finite system this is valid 
approximately, if N is large. Indeed, in this case for a uniform distribution of 
constants ipk — "01 + 27r(/c — 1)/N the computation of C\ according to (1361 
yields \C\\ = 1, arg(Q) = ipil for I = N, 2N, . . ., and C\ = otherwise, and we 
get 

7 = 1 + (1 - p- 2 )-^~ r —t-w ■ (38) 



1 - [-pe*Wi-*) 

Hence, the deviation of 7 from unity decreases with the size of the subpopu- 
lation and, therefore, can be neglected for large N. 

Using again the subpopulation index a, we summarize that for the uniform 
distribution of constants of motion and for large iV a the order parameter of a 
subpopulation is directly expressed via the WS variables: 

Z a = p a e iq>a = z a . (39) 

As a result, the first two WS equations (see, e.g., Eqs. (130}) ) become equations 
for the amplitude and phase of the local mean field. The system simplifies 
further if the forcing H is independent of ^, then the third WS equation 
becomes irrelevant. 



4-2 WS variables versus generalized order parameters 



The Kuramoto order parameter Z defined according to Eq. (jTBl or Eq. (I3TT) 
is an important quantity but it does not provide a full characterization of the 
oscillator population. Such a characterization is given by a set of generalized 
Daido order parameters [32|33|45ll3l] . A generalized Daido parameter of the 
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order m is defined according to 



iV 



or 



k=l 



2- 



w((j>)e im *d(j> 



(40) 



Clearly, Z\ is just the Kuramoto order parameter Z. The physical meaning 
of the parameters Z m is especially transparent in the thermodynamic limit: 
they are just the Fourier harmonics of the distribution of the phases and thus 
completely characterize this distribution. Using (jSJ) we obtain: 



p m e tm *>y m (p,#) = z m lm {z,a) 



where 



N 



Jm (p, vl/) = N- 1 £ 
fc=i 



p -i e i(^ fc -*)\ m 



or 



7m (p,#) 



1 + p&fa-*) 
l + p -i e ^*)\ m 



[for brevity we skip a similar expression for 7 m (z, a). 



a(i/j)dijj 



(41) 

(42) 
(43) 



It can be seen that for uniform distribution of constants of motion in the 
thermodynamic limit, i.e. for a = (27r) _1 , j m — 1 for all m. To show this, we 
again write 7 as a series and obtain 



7n 



/_ r 00 

(l + /r i e ^-*)) m i + £(_ pe ^-*) 



adip 



Computing the powers and performing multiplication we obtain a series where 
the first term is (2ir)^ 1 J^^ dip = 1 and all other terms are products of some 
functions of p and \l/ with the integrals (27r) _1 e tL ^dip = 0, where L is an 
integer. Thus, for the special case of uniformly distributed constants of motion 
we obtain j m = 1 and 

Z m = z m = Z m . (44) 



Returning back to the notations for subpopulations of oscillators, we can write 
a general relation between the local order parameters and the WS variables 

as 

Z a , m = Pae tm * a la,m(Pa, **) = C7a,m(*., <*.) . (45) 

For a particular case of uniformly distributed constants of motion this relation 
simplifies to 

Z a , m = z™ = Z™ . (46) 
We emphasize that all results of this Section are valid for the case of a pop- 
ulation with a continuous distribution of parameters as well. In this case we 
deal with the order parameters 

Z m (x,t)= r 'wl^le^ (47) 
Jo 
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and functions 7(2;, p, $). For uniformly distributed constants if) we have j m (x) = 
1 and Z m {x) = Z m {x). 



4-3 The Ott-Antonsen ansatz 



Ott and Antonsen [23] treated the Kuramoto problem in the thermodynamic 
limit N — > 00 with the help of the continuity equation ( 12T1) . Writing the 
density function w(u), <f>, t) Fourier series 



n(uj) \ 
w{u,<j),t) = < 1 + 



^/ m ( W) t)e-^ + c.c. 



.m=l 



where c.c. denotes complex conjugate, Ott and Antonsen noticed that the 
continuity equation is fulfilled if the Fourier coefficients can be expressed as 

f m (u,t) = [F(u,t)r , (48) 

where F(u,t) is some function. 

Comparing this approach with the definition of generalized order parameters 
PPl) . we see that the quantities f m are exactly these order parameters and the 
ansatz (HHi) means that 

Z m (cu,t) = [Z(cu,t)] m 

Thus, the found particular class of solutions, which we denote as the OA 
reduced manifold, exactly corresponds to the special case where the general- 
ized order parameters are expressed via the powers of the WS variable z (see 
Eqs. ( [441146]) . This holds if the distribution of the WS constants if> is uniform. 

The OA ansatz can be alternatively presented as follows. Let us consider a 
generalized order parameter Z m (u,t) of a subpopulation with the frequency 
u, see Eq. fT4"T]) . and compute its time derivative 

Jo at Jo 

here we also used Eq. ([27]). Substituting <fi = u + (He^ - H*e i(t> ) /2i we obtain 
(cf. 121): 

m 

Z m = iumZ m + —{HZ m -\ — H*Z m+ i) . 

This (infinitely dimensional) system of ODE obviously simplifies if Z m = Z m ; 
this solution exactly corresponds to the OA manifold. On the other hand, it 
corresponds to the particular solution of the WS equations for the uniform 
distribution of constants of motion 7(0;) = 1. In this case the WS equations 
and OA approach yield the same equation for the time evolution of the order 
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parameter: 

Z(u, t) = iuZ + ^(H - H*Z 2 ) . 
This important issue is illustrated in the next Sections. 



4-4 Reduced solution of a system with continuous frequency distribution 



In the very recent publication [21], Ott and Antonsen have demonstrated that 
the reduced manifold ( |48l) is the only attractive one if the globally coupled 
system has a continuous distribution of parameters. Following their ideas, we 
demonstrate how this property follows from our theory for a particular case 
of harmonic force H with frequency u; for this purpose we use the results of 
Section 13.1.41 We remind, that in the terms of WS approach, the reduced OA 
manifold corresponds to the case j(x) = 1. 

Suppose we are interested only in a global characterization of the dynamics, 
e.g., we compute the global mean field, which is obtained via a summation or 
an integration over the whole population: 

Y = J dxn(x)Z(x)= J dxn(x)-f(x)p(x)e i ' s>{x) . (49) 

Substituting here 

00 

j(x) = 1 + (1 - p(x)- 2 ) £Q(x)(-p(x)e-^))< , 

1=2 

to be compared with Eq. (1371) . we obtain 

Y = J dxn{x)^(x)p{x)e iq ' { - x " 1 = 

= f dx n{x)pe iq ' + Y.i- 1 ) 1 [ dx n(x)pe i ^C l (p l - p l - 2 )e~ m . (50) 

■* l>2 



We argue, that the sum in Eq. (1501) eventually vanishes and only the first 
term remains important. In Section r3.1.4l we have shown that for the case of 
a harmonic forcing, two types of solutions are possible for a subpopulation 
with the frequency u(x): a fully synchronous attractive solution with p = 1 
and a quasiperiodic solution described in Appendix [B] For the interval (s) of 
parameter x where p = 1, all terms in the sum in Eq. floTJ]) vanish. Thus, 
we have to consider only integrals over the intervals where p(x,t) < 1. To 
this end it is important to note that the expression under the integrals in the 
sum contains oscillating functions p, and According to (1B.9j) . the phase 
variables $, \1/ rotate with a frequency that smoothly depends on x. Hence, 
the expression under the integrals oscillates with some frequency, smoothly 
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depending on x. Thus, for large t this term rapidly oscillates in x and, hence, 
the integral over x vanishes, provided the other dependencies on parameter 
x are sufficiently smooth. Mostly important, the distribution of oscillators 
n(x) should not have singularities. Otherwise, if n(x) contains 5-functions, 
the integrals in the sum do not vanish. 

Summarizing, we expect that for t — > oo the sum in Eq. ( 150]) tends to zero and 
the global order parameter can be expressed via the WS variables p, $ only: 



Thus, the functions 7(2) become irrelevant for the macroscopic dynamics. For 
t — > 00 the mean field can be computed as if 7 = 1, for all frequencies, what 
corresponds to the OA reduced manifold. Physically, this occurs due to the 
effective "collisionless" mixing of different subpopulations that are not locked 
by the common force. This effect is similar to the Landau damping in plasmas 
and to the inhomogeneous line broadening in optics. 

We note that our argumentation treats only the case of harmonic mean field, 
because it is based on solutions of the WS equations ( jT6]fT8]) . A rigorous proof 
for arbitrary functions u and H was given recently by Ott and Antonsen [21] . 



5 Mean field coupling 

5. 1 Organization of a subpopulation: mean field coupling 

Till now we considered a general time dependent force H a , acting on elements 
of the subpopulation a, just it had to act equally on all elements. Now we 
specify this force and consider several popular models as particular examples 
of the general approach. 

Generally, the force H a can have two components. The first one arises from 
the interaction between elements of the subpopulation a itself (Fig. [p. Very 
often it is assumed that this component is computed in a mean field fashion, 
i.e. that the coupling within each pair of oscillators is the same. Hence, the 
component of the forcing due to internal interaction is proportional to the 
complex mean field (order parameter) Z a of the subpopulation (see Eq. (fl3l)). 
The proportionality factor we denote by E aa ; generally it is complex. 

The second component results from all forces external with respect to the 
subpopulation a: those are the forces from all other subpopulations, regular or 
noisy forces acting on the whole population, etc. Next usual assumption is that 




(51) 
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all elements of a equally contribute to the force acting on other subpopulations, 
and, vice versa, and therefore the effect of subpopulation b on subpopulation 
a is proportional to its complex mean field Z b and to its size n b , so that 

M 

Fa.ext 

+ E ab n b Z b . (52) 

6=1 

Here F a e3St is the sum of all other external forces, acting on a, complex con- 
stants E ab describe the coupling from subpopulation b to subpopulation a 
(the term in the sum with a=b corresponds to the first component described 
above), and n b = N b /N are relative population sizes0 At this point we em- 
phasize that coupling determined by Eqs. (1521) is only a simplest form of the 
mean field coupling which we denote as linear. The general, nonlinear mean 
field coupling will be discussed below in Section [61 and now we consider several 
examples of linearly coupled ensembles. 

First of all we derive the closed set of WS equations for interacting subpop- 
ulations with the mean field coupling. For the case of discrete set of subpop- 
ulations we complement Eqs. f l22|23|) by the expression for the force Eq. ( 152]) 
which we re- write, using Eq. ( 152"|) . as 

M 

H a = F a>ext + E ab n b 'j b z b . (53) 

6=1 

The continuous analog of this equation is 

H(x) = F ext (x) + J dyE(x,y)n(y) 1 (y)z(y) . (54) 

It complements Eqs. ( )25|26|) in case of infinite number of subpopulations. In 
these relations 7 is defined according to Eq. ( I3"3"|) or Eq. (133]) . We stress here 
that the systems ( I22|I23II53I) and (I25|26ll54p represent an exact reduction of 
the dynamical description of hierarchical populations of oscillators by virtue 
of the WS ansatz. 

As already discussed, the obtained equations significantly simplify on the Ott- 
Antonsen manifold. In Section H~2l we have demonstrated that for the uniform 
distribution of the constants of motion of a subpopulation, the corresponding 
function 7 = 1. In this case one of the WS variables decouples and we obtain 



4 For the following it is convenient not to absorb n b into the coupling constant E, 
but to keep the effect of the subpopulation size explicitely. 
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a reduced system 

^ = i"az a + \H a - Z \ui , (55) 

M 

H a = F a>ext + E ab n b z b , (56) 

6=1 

for a discrete set of interacting subpopulations, or, respectively, 

= &;(*, t)z(x, t) + \h(x, t) - ^H*(x, t) , (57) 
H(x, t) = F ext (x, t) + J dy E(x, y)n(y)z(y, t) , (58) 
for a continuous distribution of parameters. 

5.2 Example 1: The Kuramoto-Sakaguchi model 



Suppose that external forces are absent, i.e. in Eq. fl52|) F aext = for all a. 
Writing the coupling constants as E ab = e ab e ll3ab and substituting Eqs. ( IT31I521) 
into Eq. <^ we come to the model (cf. |^2|23] ): 

jjfl) i M N b 

-^r = w « + ^EE^ M4 } - A a) + Pa b ) , k = i, . . . , N a . (59) 

at iV 6=1 j=l 



If we furthermore assume that the coupling parameters are the same for all 
subpopulations, e ab = e, (3 ab = (3, then we obtain the well-known Kuramoto- 
Sakaguchi model [3Tfl6] of globally coupled oscillators: 

,(a) ^ N 



1 "a + e-YMtPi-^ +P) , k = l,...,N a , (60) 



dt N l=1 

where the summation is over the whole population containing N oscillators. 
In terms of the global mean field Y = re 10 the model reads 

Ti(a) 

~- u a + er sin(6 - 4 a) + (3) , k = l,...,N a , (61) 



The analysis of this model by virtue of WS reduction for the case of identical 
oscillators has been performed in the original WS paper [39J . 

We proceed here by discussing the Kuramoto-Sakaguchi model with a con- 
tinuous frequency distribution; its description is given by Eqs. (I25f26|54l) . For 
this model it is natural to identify the continuous variable x with frequency 
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oj. For the common effective force we obtain 



H = EY = ee il3 Y , (62) 

where Y is the global mean field (see Eq. (j49p ). Substituting the force into 
Eqs. ( }2"5|2"6"|) we obtain a closed system of equations 

dz(u,t) _ , E* 2 



iuz + -Y-—z 2 Y* , (63) 
u + lm(z*EY) . (64) 



dt ' 2 2 

da(oj, t) 
dt 



Consider now the reduced set of equations for the Kuramoto-Sakaguchi prob- 
lem. This reduced set corresponds to the case 7(0;) = 1, z(cu) = Z(cu) and 
contains only two equations (I63|62p . Next, we consider the Lorentzian distri- 
bution of natural frequencies, n(cu) = [ii{ui 2 + l)] -1 . As demonstrated by Ott 
and Antonsen [23], for this case, under an additional assumption that z(u) is 
analytic in the upper half-plane, the integral in Eq. ([4"5]) can be calculated by 
the residue of the pole at u = i; this calculation yields Y = z(i). Substituting 
this along with u = i into Eq. (16111) we obtain the OA equation for the time 
evolution of the Kuramoto mean field: 

£-(- 1 + f)r-fr*.. <«> 

This closed equation for the order parameter was first derived and solved 
in 1231; the solution is 



r(t) = R{l + 



^0 



'1/2 



\ , (66) 



where R = yl — 2/e (notice a misprint in Eq. (11) of (23])- However, gener- 
ally solutions of ( 165]) do not coincide with the solutions of the full equation 
system ( l63|64l) . We illustrate this important issue by the following numerical 
examples. 

First, we show that the solution deviates from ( )66|) if the analyticity assum- 
tion above does not hold. We perform a direct numerical simulation of the 
Kuramoto-Sakaguchi model with iV = 10 4 oscillators. The frequencies of the 
oscillators are all different and are chosen to approximate the Lorentzian dis- 
tribution. For the parameters of coupling we take (3 = 0, so that E = e is 
real. We perform two runs with the same macroscopic initial conditions for 
the ensemble, choosing Y(0) = r(0)e je< - ^ = r = 0.5, but with different initial 
distribution of phases. Practically, we introduce an auxiliary angle variable ? 
which attains iV values, labelled by index k, uniformly distributed between 
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Fig. 4. (Color online) Mean field amplitude r as a function of time, for the Kuramoto 
ensemble with the Lorentzian distribution of frequencies; the number of oscillators 
is N = 10 4 , e = 3. Curves (a) and (b) correspond to two different sets of initial 
phases, as described in the text. In the first case the evolution of the mean field 
follows theoretical solution given by Eq. (166p , while for the second case the transient 
dynamics deviates significantly from this solution. 

— 7r and 7r (end points are excluded). The frequencies of oscillators are then 
obtained as Uk = tan — and the initial phases as 



±2 arctan 



r ° tan(<? fc /2) 
1 + r 



±2 arctan 



l-r 
— U k 

1 + r 



(67) 



cf. Eq. (El); here the plus and minus signs correspond to the first and the 
second run, respectively. Using Eq. §5§ and Appendix [A] we write for the WS 
variable 

= = r„(lT^) + l±^ ( 

r (1 ±lUJ k ) + 1 =F lUk 

(Notice that p(u)k) — 1 because we have only one oscillator at each frequency.) 
Considering the obtained expression as an approximation of a continuous func- 
tion z(co), we find that the latter has a pole at w = T^i~^- Thus, the first 
case corresponds to an initial condition that is analytic in the upper half-plane, 
while the second run corresponds to initial conditions that are analytic in lower 
half-plane. The results are shown in Fig. HI we see that the transient dynamics 
of the global mean field heavily depends on the microscopic initial conditions. 
We emphasize that the result for the first set of initial conditions very well 
agrees with the solution fl6"61) . while for the second set of initial conditions the 
transient dynamics is essentially different. 

Next, we verify the validity of the theory for hierarchically organized pop- 
ulations, by simulating the same model with M = 10 4 groups of identical 
elements. All groups have the same size, i.e. N a = N g for all a = 1, . . . , M. 
Again, we always start with the same macroscopic initial conditions for the 
ensemble, choosing Y(0) = 0.5. However, the microscopic initial conditions, 
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Fig. 5. Mean field amplitude r as a function of time, for e = 3, M = 10 . In all 
simulations initial conditions have been chosen so that r(0) = 0.5. For N g = 100 
the evolution of the mean field follows Eq. (|66p . while for smaller group sizes the 
transient deviates significantly from the OA manifold. In all runs the distribution 
of the constants ip a ,k inside each group was chosen to be uniform (q = 1). 
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Fig. 6. The same as Fig. but for different distributions of constants of motion 
V> a ,fc and fixed group size N g = 200. Dotted line shows the theoretical asymptotic 
value of r. 



given by the distribution of the constants of motion ifj a ,k, differ from run to 
run. In particular, we introduce a parameter < q < 1 that quantifies devia- 
tion of the distribution of ip a ^ from a uniform one; the value q — 1 corresponds 
to the uniform distribution. In Appendix [D] we describe how one can choose 
different microscopic initial conditions while keeping the same macroscopic 
initial conditions. Here we choose the initial state in such a way, that z(oj) is 
analytic in the upper half-plane. 

First we analyze the effect of the subpopulation size N g (Fig. EJ), taking q = 1. 
Theoretically, this effect is described with the help of Eq. (138|) . The results 
confirm the theoretical prediction: with increase of N g the transient dynam- 
ics tends to the OA manifold and is nicely described by Eq. flrJrJj) . while for 
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small number of oscillator in a group, the deviations from the OA solution are 
essential. If there is one oscillator in a group, the OA solution (166]) is again 
valid because here p = 1. The effect of a non-homogeneous distribution of 
the microscopic constants ipk,a on the dynamics of the mean field, for a large 
group sizes, is illustrated in Fig. |6j One can see that the deviations from the 
OA manifold become larger as this distribution becomes less and less uniform, 
i.e. if the parameter q deviates from one. The examples of Fig. I4"t5f6l illustrate 
that although the OA ansatz yields a simple closed system of equations, these 
equations do not describe the dynamics for general initial conditions, but only 
for a special subset of them. 



5.3 Example 2: Two coupled subpopulations 



For the next example we concentrate on a model, recently studied by Abrams, 
Mirollo, Strogatz and Wiley [IS] . They considered two identical subpopula- 
tions of the same size, i.e. Ui = u 2 = u (without loss of generality we set it to 
zero) and Aq = N 2 , but the coupling within a subgroup differs from the cou- 
pling between the subgroups: Su = e 22 = £12 = £21 — v 7^ an d fi a b = ft- 
The equations are: 



d4 a) 



lo + Im 



(jiZ a + uZ b )e^-*^ 



dt 

where a = 1,2. The WS system (I22II23I) for this setup reads 



(69) 



M4Hi) , (71) 



- z h • o» 

dat\ 
~dt 

^ = lm(zlH 2 ) , (73) 
#i,2 = (^i, 2 + vZ 2A y? , (74) 

and the relation between Zx j2 and is given by Eq. f )32|) . By applying the 
OA ansatz, i.e. by setting Z 12 = 24,2 we obtain a set of equations, originally 
derived in Ref. [18J. Analyzing these equations, Abrams et al. have obtained 
an interesting solution where one subpopulation is completely synchronized, 
\zi I = 1, while the other one is only partially synchronized, \z 2 \ < 1. Moreover, 
this partially synchronous state can be either steady, z 2 = const, or time- 
periodic, i.e. z 2 is a periodic function of time. These regimes are called chimera 
states. 
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Fig. 7. (Color online) Simulation of ensemble (|69[) for N = 64, (3 = tt/2 — 0.1, 
/x = 0.6, v = 1 — [j, = 0.4, and different distributions of the microscopic constants 
ip k , controlled by the parameter q (see Appendix [D]). Note that the distribution of 
constants Vu,fc is irrelevant since the first subpopulation is completely synchronized. 
The case q = 1 (marked by plus) corresponds to the OA manifold, here the mean 
field is constant. For q = 0.9, q = 0.7, and q = 0.5 one observes time-periodic states 
represented by limit cycles in the complex plane Z2 (red bold, green solid, and blue 
dotted curves, respectively). 

The model of Abrams et al. serves as a good illustration of the usefulness of the 
above described approach based on the exact WS theory. A complete descrip- 
tion of the dynamics for arbitrary initial conditions is given not by the OA 
equations, but by system (I70H74||321) . Correspondingly, the additional equa- 
tions generally lead to an additional time-periodicity for chimera states |40j : 
a steady-state solution becomes time-periodic (Fig. [7]), and a time-periodic 
state becomes quasiperiodic (Fig. El). We notice, that in this case the solu- 
tions do not evolve towards the OA manifolds, because the distribution of the 
oscillators' parameters is not continuous. 



6 Nonlinear mean field coupling 

6.1 General nonlinear coupling 

Till this moment we restricted our consideration to the case when the force, 
acting on the oscillator subpopulation a, is a linear combination of the local 
mean fields Z h of all other subpopulations b (see Eq. (152"]) ). We denoted this 
coupling as linear. Generally, this force can depend on the Daido generalized 
complex order parameters |32,34j, see Eq. (140)) . For a general nonlinearity we 
can expect that H a contains arbitrary combinations of Z^ m . However, in a 
physically reasonable model some restrictions appear. 
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Fig. 8. The same as in Fig. but for [i = 0.65 and v = 1 — n = 0.35. Now on the 
OA manifold, i.e. for q = 1, the dynamics is periodic (a), while for more general 
initial conditions parameterized by q = 0.9 (b), q = 0.7 (c), and q = 0.5 (d) the 
dynamics is quasiperiodic. 

To discuss this important issue let us first consider an isolated subpopulation, 
say a, and recall that derivation of the phase models of type ([T]l3]) incorporates 
an averaging over oscillation period (see, e.g., [H] ). It means that the mean 
field forcing H a can include only the terms having the frequency w u, i.e. the 
terms like 

Za,l i 2/ am Z a.l—m i ' 'a,m^ a,l^ a,l— m— 2 j • • ■ ; 

where the sum of lower indices is one (we remind that Z_ m = Z^). Let us now 
consider a particular case when the nonlinear mean field coupling is determined 
by the first order parameter Z a = Z a ^ = r a e lBa only. Remarkably, this case 
corresponds to OA reduced manifold since the latter implies Z m = Z m . The 
forcing then takes the form: 

H a = h,Z a + h 3 \Z a \ 2 Z a + h 5 \Z a \ 4 Z a + ... . (75) 

Denoting \h x \ = e and 1 + ^\Z a \ 2 + ^\Z a \ 4 + . . . = A(r a , e)e i/3(ra ' e) we obtain 

H a = eA{r a) e)e^ ir ^Z a (76) 

and 

= u a + eA(r a ,e)r sin (Q a - <p k + 0(r a ,e)) . (77) 

This nonlinear model was suggested and analyzed in [2TH22] ; its interesting 
dynamics is discussed below. 
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In case when the subpopulation a interacts with the other populations and 
external fields, Eq. (!75|) includes additional terms. The terms, describing inter- 
action with other subpopulations can contain any combination of -Z& jm , with 
b = 1, . . . , M, however, only the terms having the frequency ~ u are resonant 
and therefore essential for the dynamics. This is the consequence of the fact 
that our basic model ([1]) is not general but contains only the first harmonics. 

6.2 A minimal model of nonlinearly coupled oscillators 

Our generalization [2Tf22] of the model ( |60l) accounts for a possible nonlinear 
response of the oscillator to the forcing. It means that the effect of the large 
force is not just an "up-scaled" effect of the small one, but can be qualitatively 
different. For a more detailed explanation of this concept, let us consider one 
oscillator, influenced by a harmonic force with the amplitude 6 and phase 
and let the interaction be described by the sine-function, so that the equation 
for the oscillator's phase reads: 



Parameters A and (3 determine the response to the forcing. So, e.g., if A cos (3 > 
0, the force with a frequency close to u results in a stable in-phase synchro- 
nization of the oscillator. Consider now an ensemble of identical oscillators, 
coupled via the mean field. Comparing Eq. f!78|) with Eq. flBTJ) we identify S 
with er and O with the phase of the mean field. If the oscillators are syn- 
chronized, the mean field has the same frequency and the phase as each of 
them, and if the above condition A cos (3 > is fulfilled, then this synchronous 
regime is stable. Otherwise, if Acos [3 < 0, the in-phase synchrony is unstable 
and the oscillators remain asynchronous, i.e. they have different phases (the 
frequencies are the same since the oscillators are identical). 

Nonlinearity of the coupling means that the parameters in Eq. f!78|) can depend 
on the amplitude of the force, A = A(5), (3 = /3(<5)H] In this case we can 
expect interesting dynamics to occur. Suppose the factor A cos (3 is positive 
for small 5 but becomes negative when 5 increases and achieves some critical 
value 5 C . Then the synchronous state becomes unstable and the oscillators 
tend to desynchronize. However, this immediately reduces the mean field, i.e. 
the amplitude of the forcing 5, and the synchronous state becomes stable 
again, the oscillators tend to synchronize, what increases the mean field. As 
a result of the counter-play of these two tendencies, the system settles at the 
border of stability, exhibiting partially synchronous dynamics: the (identical) 
oscillators are not synchronized, because synchrony is unstable, but they are 

5 Generally, functions A and (3 can depend not only on the product of r and e, 
but on both variables, i.e. A = A(r, e), j3 = (3(r,e), see [22] for details. 



= u + A5 sin(6 - fc + (3) . 
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also not completely asynchronous, because this state is unstable as well. As a 
result, the oscillators distribute around the unit circle, forming a bunch. Two 
different dynamical states appear depending on whether A or cos /3 becomes 
negative with an increase of 5. In the former case this bunch remains static 
in the frame, rotating with the oscillator frequency u osc . In the latter case, 
this bunch rotates with respect to the mean field, i.e. the oscillators have 
frequency that is generally incommensurate with the frequency of the mean 
field. Furthermore, the bunch can also "breath" so that the mean field is 
modulated, see [2"1~P2] for more details. 

Let us now take an infinitely large population with a continuous frequency 
distribution. Consider a particular case of the homogeneous coupling when all 
oscillators are driven by the same force. Then we should omit the index a in 
Eq. f!75|) and substitute there Z a by the global mean field 

Y = re i0 = J n{uo)Z{uj)dLo , H = eA(r, e)e^ {r ' £) Y (79) 

This corresponds to the following microscopic equations 

<p k = u k + eA(r, e)r sin (6 - <p k + /3(r, e)) . (80) 

Next, we consider the infinitely large system and want to write the corre- 
sponding WS equations. We emphasize that in the derivation of Eqs. (j63|64p 
we did not assume that E = const; we only used the assumption that the 
coupling is homogeneous, i.e. that E is independent of the frequency. Hence, 
we just have to substitute in Eqs. (I63II64P E with eA(r, e)e l ^ r,e \ This yields, 
together with Eq. (179]) . a system of integral equations with nonlinear coupling. 
This full system is still rather difficult to analyse, therefore below we perform, 
like in Section 15.21 a further simplification allowing us to obtain an analog of 
Eq. ([65]). 



6.3 Nonlinearly coupled ensemble with the Lorentzian frequency distribution: 
Theory 

In this Section we exploit the general theory to analyze in details the dynam- 
ics of ensembles with global nonlinear coupling (see Eqs. (180]) ) and Lorentzian 
distribution of frequencies. Looking for the asymptotic solutions we follow 
the argumentation of Section H] and consider the reduced dynamics, corre- 
sponding to the case 7(0;) = 1 and, respectively, to the OA reduced manifold. 
Furthermore, following |23j we consider the Lorentzian distribution of natural 
frequencies, n(u) = [vr(u; 2 + l)] -1 , and assume that the field z(co) is analytic 
in the upper half-plane. Then, in a way similar to the derivation of Eq. (I65p . 
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we obtain 

dY = ( eA(r,e)e^e) \ _ £ A(r,e)e^) 
dt \ 2 J 2 

Below we verify the validity of this equation by numerics; in particular we 
confirm that the asymptotic solutions are confined to the OA manifold. How- 
ever, the basins of attraction of these solutions depend on the distribution of 
the microscopic constants of motion. 

Separating the real and imaginary parts of Eq. flHTl) we obtain equations for 
the amplitude r and frequency Q of the mean field: 

dv eA 

— = -r + — r(l - r 2 ) cos , (82) 
^ = fi = ^(l + r 2 )sin/3. (83) 

Stability of the asynchronous state r = of the ensemble is determined by the 
condition 



df 
dr 



= -1 + | A\ r=Q cos P\ r=0 = -1 + | A cos p < , (84) 

r=0 - " 



what yields the value of the critical coupling 

2 



A cos j3 Q 



J5) 



For e > e cr the ensemble exhibits a synchronous state with the mean field 
amplitude < r < 1. @ The latter is determined from the condition r = 0, 
which yields the equation 

e(l -r 2 )A(r) cos P(r) = 2 . (86) 

If the solution of this equation for particular functions A, P is found (most 
likely, numerically), then Eq. (183]) provides the frequency of the mean field Q. 
Note that the collective oscillation arises with the frequency Q \ £ =e cr = tan(/3o). 
Below, we illustrate the theory by two particular choices of the amplitude 
A(r,e) and phase P(r,e) dependencies. 



6 It can be shown that for a physically reasonable model (see Section I6.ip . the 
order parameter grows close to the transition point as r ~ y/e — e cr , as in case of 
the standard Kuramoto model. 
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Fig. 9. Mean field amplitude r as a function of coupling strength e for the nonlinear ly 
coupled ensemble with A = 2 — er 2 and j3 = 0. Black bold line: theory. Symbols 
show simulations of hierarchical ensembles with M = 1000 groups of 20 oscillators, 
for different initial conditions (see text). Red line corresponds to simulation of a 
population, where all elements have different frequencies, without groups, and with 
a uniform distribution of (j). 

6.4 Nonlinearly coupled ensemble with the Lorentzian frequency distribution: 
amplitude nonlinearity 

First we consider the impact of the amplitude function A, by setting A(r) = 
fi — er 2 , [5 = 0. Equations ( )83l) and ( 1851) yield Q = and e cr = 2//x. From 
Eq. (1861) we obtain a quadratic equation for r 2 



It is easy to see that for e > we have /(±oo) = 00, f(0) = sfi — 2, and 
f(l) = —2. Hence, for fi > 0, the solution in the [0, 1] interval exists for 
e > e cr = 2//i, it is given by 



The dependence of r on e for \i — 2, e cr — 1 is shown by bold line in Fig. |9j 

We verify the theory by simulation of ensemble dynamics for M = 1000 groups 
of N a = 20 oscillators. For this goal, we prepare different initial conditions, 
corresponding to uniform and nonuniform distribution of ipk,a] these distribu- 
tions are parameterized by parameter q, so that q = 1 and q < 1 correspond to 
uniform and non-uniform distributions, respectively (see Appendix [D]) . Next, 
we simulate the ensemble of iV = 10000 oscillators with different frequencies 
(in other words, each group contains only one oscillator); in this case we take 
for initial conditions a nearly uniform distribution of phases <pk- We see, that 
the results, shown in Fig. [9J demonstrate a good correspondence to the theory. 



e 2 r 4 - e(e + [i)r 2 + e/x - 2 = f(r 2 ) = . 
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6.5 Nonlinearly coupled ensemble with the Lorentzian frequency distribution: 
phase nonlinearity 

Now we analyze the effect of the dependence (3 = (3(r), by setting A — 1, 
/5 = /5o + £ 2 r 2 . For the chosen particular function we have e cr = 2/ cos (3q and 
Eq. ( iSljj) yields an equation for r 2 : 

e(l - r 2 ) cos(/3 + e 2 r 2 ) - 2 = /(r 2 ) = . (89) 

It is easy to see that /(0) = (e — £ cr ) cos/?o > and /(l) = —2, hence there 
always exist at least one solution. (We remind that (3o = const, |/?o| < tt/2.) 
Numerical analysis of Eq. (189]) shows that the number of its roots increases 
with e. Thus, the system exhibits multistabilty. The corresponding bifurcation 
diagram in the parameter plane e is shown in Fig. [TTJ1 Dependencies of the 



as/3s 




A) 

Fig. 10. Multistability in the nonlinearly coupled ensemble with A = 1, /3 = fto+e 2 r 2 . 
Red dashed line shows critical coupling e cr = 2/ cos/?o; inside the domain, deter- 
mined by this curve, the asynchronous state is unstable. Labels as, s, and ns mean 
asynchrony (the state with r = is stable), synchrony (r > 0), and coexistence of 
n synchronous states, respectively. Label as/ns means coexistence of asynchronous 
and n synchronous solutions. 

mean field amplitude and frequency on the coupling strength for (3 = are 
shown in Fig. [TTJ 

To verify the theory we again perform a direct numerical simulation of an 
ensemble with M = 5000 subpopulations of iV a = 20 oscillators each, for 
(3q = 0. The numerical results are shown by symbols in Fig. [TTJ 

Finally, we demonstrate that although the stationary dynamics of the system 
corresponds to uniform distribution of microscopic constants ip (i.e., to q = 
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Fig. 11. Illustration of the multistability in the nonlinear ly coupled ensemble with 
A = 1, f3 = (3q + e 2 r 2 , for (3q = 0. Three branches of the solution of Eq. (|89p for the 
mean field amplitude are shown by different colors in (a). Corresponding solutions 
for the frequency of the mean field are shown by the same colors in (b). Numerical 
results (see text for details) are shown by symbols. 

1), the transient dynamics does depend on the distribution. In other words, 
the attractors of the multistable system can be obtained by the simplified 
theory, see [23] and discussion above, but their basins of attraction depend 
on the distributions of ip. In numerical experiments, we simulate an ensemble 
containing M = 5000 subpopulations of 20 oscillator each, for e = 4.5, taking 
Po = 0.52 and different values of q. The results shown in Fig. [T2l demonstrate 
that starting from the same macroscopic initial conditions, the system can 
evolve to different attractors, depending on the microscopic constants. 



7 Conclusions and outlook 

The main goal of this paper was to provide a generalization of the powerful 
Watanabe-Strogatz theory on the heterogeneous populations of phase oscil- 
lators. We have formulated the Watanabe-Strogatz equations for a general 
hierarchically organized ensemble, and have examined limiting cases of infi- 
nite populations. Remarkably, there exist two possible thermodynamic limits: 
in the first one we treat a finite number of infinitely large populations, whereas 
in the second case we deal with a system with a continuous distribution of 
parameters, e.g., of frequencies. The derived equations provide an exact reduc- 
tion of the dynamics; in many cases the problem under consideration becomes 
low-dimensional. We have analyzed the derived equation in several important 
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Fig. 12. Evolution of the mean field for same macroscopic initial conditions 
(r = 0.52, = 0, shown by a cross) but for different distributions of the con- 
stants of motions, parameterized by q (see text). Note that both attractors of the 
system (limit cycles with the radius ~ 0.23 and ~ 0.6) correspond to the theory, 
developed under assumption of the uniform distribution of the constants of motion 
(cf. Fig. [TT1) . However, transient and basins of attraction depend on the distribution 
of the constants of motion. 

cases, including the Kuramoto-Sakaguchi model and the model of two coupled 
populations with chimera. Noteworthy, the reduced equations are valid both 
for linear and nonlinear coupling; in the latter case the approach has allowed 
us to describe multistable synchronous dynamics. 

Next, we have thoroughly studied a relation between the Watanabe-Strogatz 
theory and the recent Ott-Antonsen ansatz and have demonstrated that the 
latter corresponds to a particular choice of initial conditions for the ensemble. 
To be exact, the OA approach corresponds to the case when the constants of 
motion in the Watanabe-Strogatz ansatz are uniformly distributed. Although 
the Ott-Antonsen equations are much simpler than the full Watanabe-Strogatz 
system, several examples considered have shown that they provide only asymp- 
totic solutions, whereas the transient dynamics and the basins of attraction of 
these solutions depend on the choice of initial conditions. (See jl?] for another 
example of nontrivial transient dynamics off the OA manifold.) 

Finally, we would like to mention that the approach presented opens new 
perspectives in analysis of such long-standing problems as finite-size effects and 
the effects of a common external noise on oscillator ensembles. Also application 
of it to systems with delayed coupling appears promising. 

We thank A. Politi and E. Ott for useful discussions, and E. Ott and S. 
Strogatz for communicating their works prior to publication. The work was 
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A Watanabe-Strogatz equations in new notations 



According to Watanabe and Strogatz [38|39] . the system of A^ > 3 identical 
globally coupled phase oscillators 

<U>1 ' w(t)+A(t)sin(f(t) -<t> k ) = u{t) + g{t) cos <p k + h{t) sin <p k (A.l) 



dt 



admits a low- dimensional description. For arbitrary functions of time uj(t), 
g(t), and h(t), this A^-dimensional system is completely described by three 
global variables plus constants of motion ip k , k — 1, . . . , N, which obey three 
additional constraints, so that A^ — 3 of them are independent. The "global 
phases" \I/ and $ and the global "amplitude" < p < 1 obey the WS equations 

p = -(1 - p 2 )(gsm$ - /icos$) , (A. 2) 



p-ty = -Jl - p 2 (gcos§ + /isin$) , (A.3) 

p<3> = —g cos $ - h sin <J . (A.4) 

The solution of the original system (lA.ip can be recovered via the following 
transformation: 



We perform the variable substitution p, $ — » p, $ according to 

p = 2P - , # = ^r + 7T, $ = $ + 7T. (A. 6) 

1 + p 2 



Rewriting the r.h.s. of Eq. (lA.lj) as u(t)+lm \Z(t)e l( ^ k j with obvious relations 
Re(Z) = —h(t), lm(Z) = git), we obtain the system of WS equations (J3E1) in 
new variables. The transformation ( 1A.5I) now takes the form 

, ( <j> k -® \ I- p. ( *pk-* \ (A7 , 
tan {—^) = —p tan {—^) ■ (A ' 7) 

It is convenient to re-write this transformation in the exponential form, using 
the following identity: 

fa _ l + itan(a/2) _ [1 + itan(a/2)] 2 



l-itan(ar/2) l + tan 2 (a/2) 

,a / 2 a a\ 

cos — • 1 — tan — h 2i tan — = cos a + % sin a 
2 V 2 2 
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With the help of this identity we write: 





* 


2 






2 




V'fc- 


-* 


2 







2 1+p 

(1 + p) cos -p) sin 

~ (1 + p) cos + i(p - l)sin 

per i{ip k -^)/2 + e i(V« fc -*)/2 
~~ p e i(Vfe-*)/2 _|_ e -i(i/. fe -*)/2 ' 

what yields the desired transformation (JSJ). 

B Dynamics of WS variables in an external field 



Here we present the exact solution of Eqs. fll9H2ip for the case \u> — v\ > Hq. 
The system has one steady state A = and 



-(u - u)^ 1 + J(uj - u) 2 H 2 - 1 for (u-u) > H , / s 

Po = { (B.l) 

-(w - z/)^ 1 - t/(w - ^) 2 # - 1 for {u-v)< -H . 



Equation (I2T]1 then yields 

d = w — v — H p = 2(uj — v) =F v/(w — z/) 2 — i/g • 
Hence, \I/ rotates with the frequency 



K = 3v-2u±yf(GJ-v) 2 - H$ . (B.2) 

Consider now solution with p ^ po and A ^ Ao. Introducing new variables 
u = pcosA, v = psinA, we write the first two equations as 



ii = — [OJ — v)u — HqUV , 

• f ^ -u^ri-u 2 2, ( B -3) 

w = [uj — v)u + ~2~(1 + u — v ) i 
with v o = and mq = Po- Introducing w = u — mq we rewrite (1B.3[) as 



w = — (tu — v + H uq)v — H wv 



■ ( xff ! ^0, 2 2x ( R4 ) 

f = [OJ — v + H u )w + ~~^~{ w — v ) ■ 

Using an ansatz p = v(v 2 + w 2 )" 1 , q = w(v 2 + w 2 ) -1 and denoting 



K = u - v + -f^o^o = Ty - z/) 2 - i?o (B.5) 
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we simplify the equations to 



P = Kg+ T' (B.6) 
q = —up . 



The solution of the latter system is 



p = A cos(k(£ — to)) , 

H ( B - 7 ) 
q = Asm[K{t - t )) - — . 

2k 



Transforming back to u, v we obtain 



.i Mill hi / — I . i I I — 

u = p COS A = Xq + 



Asin(ft(t- t )) 



j\2 i £i _ AHp 

4k 2 k sin(re(t— to)) 



Acosfftft — tn)) 

v = p sin A 



^2 i £o. _ 

4k 2 /tsin(re(t— to)) 



In these variables Eq. ( |2T|) takes the form 



(B.8) 



d = cu — v — H u (B.9) 



and is readily solved by substitution of (1B.8|) and integration. As a result we 



obtain a quasiperiodic solution: variables p and A oscillate with the frequency 
k (it means that $ has the frequency k + v) and \l/ = $ — a possesses an 
additional frequency, found via integration of (IB. 91) . 
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C Variable transformation for continuity equation 



We perform transformation of variables in Eq. ( 127)) . using Eq. 



dw d . . dw dr dw dip d . . dip 



d dip 
dr \ d(p 



d ( dip\ dip 
Ihp \ a d0 J ' ~dt 



d ( dip\ ( dip\ dv 
dip \ d<p I \ dcp I dip 



dip 



da dip d f dip \ 
dr dcp dr \dcp J 



da dip d fdip\ 
dip dcp dip \dcp J 



dip 

~dt 



+ 



da dip d I dip" 
dip dcp dip \ dcp 



v + a 



dip dv 1 dip 
dcp dip J dcp 



da dip J d I dip\ d ( dip\ ( dip dip\ ( dip\ dv 
drdcP +a | Ih [dcp J + Ihp Vdcp J [~dt + V ~d4> J + [ !hp ) d^p 



da j dip ( dip dip s 
dip \ dcp \ dt dcp 



(C.l) 



Let us demonstrate that the coefficients at a and — — vanish if p, $, and \I/ obey 

dip 

dip dip 

the WS equations. For this goal we first compute — — h v-—. It is convenient 

dt dcp 

to use the notations / = e l ^~®\ c = e l ^~^\ Resolving Eq. ()5j) with respect 
to ip, we obtain 

ip = V - i\n(f - p) + iln(l - pf) . (C.2) 
Taking the derivative and re-arranging the terms, we obtain 



dt m J (f- P )(l-fp) 



$ + i 



(f-p)(l-fp)> 



(C.3) 



Using e 



e ^ j j = e we obtain in new variables: 



v = lo + Im 



He- 49 /* 



(C.4) 



Next, from Eq. (1C.2I) we compute, using 



df 



if- 



?tu\ (l~P 2 )/ 
dcp [(p> (/-p)(l-p/) 



(C.5) 
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Substituting into Eq. (1C.3|) the derivatives via the r.h.s. of the WS equations 
and using Eqs. (lC.4IIC.5j) . we obtain after tedious but straightforward algebra 



dip dip 

-oT + V d6 = °- 



(C.6) 



Hence the coefficient at — — = and the coefficient at a reduces to 

dip 



d ( dip \ (dip\ 2 dv 
^\d4>) + \d4>) lhp =Q 

To compute Q, we first substitute in Eq. (lC.5j) / 
obtain, after straightforward manipulations, 



p + c 
pc+1 



from Eq. (jSJ) and 



diP = {p + c){p + c*) = pc + pc* + 2 _ | 



1-p 2 

Derivation with respect to time yields 
d ( dip\ d ( dip\ ip(c* — c 



dr 



dt\d(p 



1-p 2 



1-p 2 



■ (l + p 2 )(c + c*)+Ap 



(1-P 2 ) 2 



dc ■ dc* 

Here we used — = —icty, = ic*fy. Next, we compute 
dt dt 



dv 
dip 



Im 



He 



-i<s> d pc+1 
dip p + c 



(p 2 - l)Re 



cHe 



(p + c) 



(C.7) 



(C. 



(C.9) 



Using the obtained expressions (1C.7HC.9I) . we show, after tedious but straight- 
forward manipulations, that Q = if ^ and p obey the WS equations. 



Thus, we demonstrate that the r.h.s. of the continuity equation Eq. ( 1C.1I) 

simplifies to — — — — and is therefore valid if <r(cu, ip) is a stationary distribution. 

dr dd) 



D Choice of initial conditions for simulation of hierarchical popu- 
lations 



Our goal is to choose different microscopic initial conditions, i.e. initial values 
for oscillator phases, but keep the same macroscopic initial conditions, i.e. the 
amplitude of the mean field. For this goal we proceed as follows. For each 
subpopulation with the frequency u a we take ip a uniformly distributed along 
the arcs [(1 - g)f , (1 + q)\ ] and [-(1 + q)\ ], -(1 - g)f ], as shown in Fig. EH 
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Here < q < 1 is a parameter quantifying deviation of the distribution from 
a uniform one; q = 1 corresponds do a uniform distribution, with q — > the 
distribution collapses to two points. For this construction, the subpopulation 
N a should be an even number. Note that this choice of ip a ^ satisfies constraints 
f|T2l and S'/'a,!; = 0. The initial values of the oscillator phases a ,fc(O) are 
obtained from ip a ^ according to Eq. flA.71) . 




Fig. D.l. Illustration of the special choice of the constants of motion ip^, here for 

q = 0.8 and N a = 10. The points are distributed along two arcs of length qir each; 

7r qir 
angle a = -(1 - q), angle P = —. 



Now we show that with a special choice of the initial values of the WS variables 
we can ensure the same initial value of the mean field, independently of the 
parameter q. These special values are $ a = 0, p a = p , and ty a = ||a. In order 
to compute the initial value of the Kuramoto mean field Y(0) = roe 1 ® we 
write the discrete version of Eq. fj49|) for t = 0: 

M 

r o e* 00 = p n ^la 

0=1 

/ M 

= po ( n » + i 1 

\a=l 

= Po , 

with account that all groups are of equal size, n a = n. Thus, taking different 
values of the parameter q and fixing other parameters we obtain the same 
macroscopic initial conditions (i.e. for the mean field), whereas the initial 
conditions for individual oscillators are different. 



M 



Po 2 )E^E^(-Po)'e" 1 



2nal 
M 



a=l 1=2 
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